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A Technique for Accelerating Iterative Convergence in Numerical 
Integration, with Application in Transonic Aerodynamics 


E. DALE MARTIN 

NASA Ames Research Center 

Moffett Field, California, USA 


Summary 

A technique is described for the efficient numcrica) solution of nonlinear partial differential equa- 
tions by rapid iteration. In particular, a special approach is described for applying the Aitken acceleration 
fomiula (a simple Fade approximant) for accelerating the iterative convergence. The method finds the 
most appropriate successive approximations, which are in a most nearly geometric sequence, for use in 
the Aitken formula. Simple examples are given to illustrate the use of the method. The method is then 
applied to the mixed elliptic-hyperbolic problem of steady, inviscid, transonic flow over an airfoil in a 
subsonic free stream. 

1. Introduction 

The numerical solutions of nonlinear partial differentia! equations such as those governing fluid 
flows frequently are obtained most efficiently by iterative methods. The rate of iterative convergence of 
the method chosen is an important consideration, and various means of accelerating the iterative conver- 
gence have been useful. 

One popular device for accelerating convergence of a sequence of numbers such as provided by 
iteration is Aitken’s extrapolation formula (or A- process) [ 1 ] , whose use is described in most books on 
numerical methods [2 ] and which is identified P-.S} as a simple Fade approximant if the successive 
iterates are partial sums of a power series. Shanks (3 ) provided generalizations of Aitken’s transformation 
and studied their use. In [6] Wynn gave a simple algorithm for rapid computation of one of the non- 
linear transforms studied by Shanks, and later Wynn [7J discussed application of this acceleration tech- 
nique to vector and matrix problems, including application to boundary-value and initial-value problems. 

The present paper describes a special technique for applying the Aitken extrapolation fonnula for 
accelerating iterative convergence in the numerical solution of partial differential equations. The method 
was first introduced and used in fSJ and then used in a modified form in [9J with additional results 
given in [9,10] . Although the application to be discussed is in a numerical finite-difference solution, the 
general method applies equally well, for example, to analytical solutions or to numerical solutions by 
finite-element methods. The use of the simple Aitken formula with three successive iterates is emphasized 



(even though the elegant e-algorithm of Wynn with longer sequences could be used), because the eventual 
applications are expected to be those numerical problems requiring significant computer storage. The 
Aitken formula, using only three iterates, requires less storage than other forms of the e-algorithm. 

Often the use of the Aitken formula with iterates obtained arbitrarily by successive approximations 
does not lead to a significantly improved approximation. However, because Shanks [3] showed that the 
formula works best if the sequence is “nearly geometric,” the present approach seeks to obtain successive 
iterates that are in a nearly geometric sequence. (Because of the work of Shanks in popularizing the 
Aitken formula and his valuable demonstration of the special applicability to “nearly geometric sequences,” 
our past work has referred to the simple extrapolation formula as the “Aitken/Shanks formula.”) The 
sequence of approximations can be most nearly geometric if obtained from a power-series construction. 
Therefore, the basis of the present approach is the construction of successive approximations derived 
from formal power-series expansions to obtain as closely as possible a nearly geometric sequence. The 
technique is based on the concepts of perturbation-series expansions (in the sense of Poincare; see 
Bellman [11]). An artificial parameter is introduced in such a way as to obtain three problems to solve 
for terms of a nearly geometric series, for use in the Aitken/Shanks formula. Expansion in powers of an 
artificial parameter has also been considered by Genz [5] to develop a mathematical proof (unknown by 
the authors of [8j at that writing), but the central idea in the present approach is that the artificial- 
parameter expansions are used, in combination with an “artificially extended f^orm” of the equations to 
be solved, as a device to determine most appropriate successive approximations. This technique produces 
the nearly geometric sequence of solutions, even in nonlinear problems. The previous application of the 
Aitken/Shanks transformations to acceleration of iterations in numerical integration by Wynn [7] used 
simple straightforward iterations. The results of such a procedure with use of only the simplest accelera- 
tion formula are described below for an example problem and are compared with the present method. 

The present approach based on perturbation series requires that complete perturbation solutions be 
available on the entire computation field (or entire domain of the equations) at each iteration. This con- 
cept therefore adapts well to a finite-difference method using “direct elliptic solvers” [12—15] in the 
iterative procedure to determine the solution simultaneously at all points on the entire computation field 
(rather than in successive traverses over the field as in a point- or line-relaxation method). Such methods 
have been referred to as “semidirect” [8-10]. 

After several simple examples to illustrate the method, it is applied to the problem of inviscid flow 
over an airfoil in a subsonic free stream, including conditions for which the flow equations are of mixed 
type (elliptic in an outer region, with an embedded hyperbolic region and a shock wave). This transonic- 
aerodynamic-flow problem has also been treated by Hafez and Cheng [16] using the Aitken/Shanks 
acceleration formula, but in a quite different way, in combination with a line-relaxation method. 



2. General Formulation of Method 

Consider the general partial-differential or difference equation system and the accompanying boun- 
dary conditions represented by 


LU - F(x) = NU in R, (2.1) 

BU = G(x) on B, (2.2) 

where U = U(x) is a vector function of the position vector x, L is a separable, linear, elliptic differen- 
tial or difference operator, F(x) is a given vector function and N is a possibly nonlinear operator such 
that the operation NU is a vector of the same dimension as IJ and has components that may involve U, 
X, and derivatives of the components of U with respect to the components of x. Assume for simplicity 
that B is a linear operator The boundary condition (2.2) is applied on B, which includes all appropri- 
ate boundary segments of the domain R. l-or illustration of this notation and of the method, simple 
one-dimensional examples are given in the next section. Examples treated in the earlier version of [8] 
included (i) the scalar Laplacian as L with a scalar, i//, as U, and (ii) a Cauchy-Riemann operator 
matrix as L with two components of U. denoted as u and v. The right side of (2. 1) can be compli- 
cated and can make the ec|uation system hyperbolic or parabolic in some regions [8-10] . 

In the formulation of a problem to be solved, L and F(x) are chosen judiciously and may be the 
result of “scaling and shifting" transformations [ 1 7,9] for increasing the rate of iterative convergence or 
of addition of terms f9,10| for stabilizing iterations. For treatment with additional terms, an extended 
Cauchy-Riemann solver for use in present calculations has been described in [ 18]. 

In the methods to be discussed for the iterative solution of eqs. (2.1) and (2.2), suppose U j(x), 
U 2 (x),, and U 3 (x) are successive approximations to U(x) in R Let u(x) and Up|(x) be respectively 
each a single scalar component of the vectors U(x) and Upfx) (n= 1,2,3). Then one form of the 
Aitken/Shanks extrapolation formula f 1 ,3 1 for an improved approximation ti*(x) to u(x) is 


u*(x) 


U IU3 - UT“ 
u I - 2 u2 + 113 


(2.3) 


Application of the formula in this way to individual components of U at each x separately is referred 
to by Wynn [7] as use of a “primitive inverse” of the e-algorithm. Wynn concludes that use of the 
primitive inverse is competitive with use of other more complicated inverses. The work of Hafez and 
Cheng ( 16] considers coupling of the matrix elements in the numerical solution, which is related to the 
more complex inverses of the e-algorithm. 

2.1 Artificially extended equation. For obtaining power-series solutions to (2.1) and (2.2) that are 
most appropriate for use in the Aitken/Shanks extrapolation formula, it has been found convenient to 
artificially extend eq. (2. 1) by inserting both an artificial parameter e and an “initial approximation,” 
Uq(x), to U(x) as follows. Let 


LU-F(x) = (l-e)NUo + eNU in R 


(2.4) 


- 7 - 


along with condition (2.2). Note that the solution U to (2.4) with (2.2) depends on e (as well as on 
the specified function Uo(x) ); U = U(x,e). However, at e=l. the solution to (2.4) with (2.2) is the 
same as the solution to the original equations (2.1) with (2.2). Furthermore, if Uo(x) is close to the 
solution U(x), then (2.4) is nearly the same as (2.1) and the solutions then are nearly the same. Thus, 
either of the conditions e=l or - Li makes (2.4) the same as (2. 1 ). Both of these facts can be used 
to advantage in the methods to be discussed. 

2.2 Method 1. The simplest iteration scheme is a straightforward method of successive appro.xima- 
tions. Althougli this method can be combined with use of a relaxation parameter (see [8,9] ), for simpli- 
city here we omit that useful device. If we let e = 0 in (2.4) and define Uo(x) as a previous iteration, 
we obtain the following equations for the iterative solution denoted as Method 1(a): 

LU„-F - NUn_| in R, (2.5a) 

BUn = G(x) on B, (2.5b) 

where subscript n denotes iteration number. 

If, as is frequently done, the Aitken/Shanks formula is used to attempt to accelerate the convergence 
of the iteration, we denote as Method 1(b) the solution of (2.5) for three successive iterates and substitu- 
tion of the results for one component of each U,j into (2.3). (This designation of Method 1(b) is useful 
for a comparison in an example problem below.) 

2.3 Method 2. The new- approach for applying the Aitken/Shanks formula, first introduced and 
used in [8] and in a modified form in (9] , is referred to as Method 2. The two versions are called, respec- 
tively, Methods 2(a) and 2(b) for later convenience. 

Consider the solution to (2.4) with condition (2.2). The solution evaluated at e = 1 is a solution to 
(2. 1 ) with (2.2). The specified Uq(x) can be used as an initial approximation to U. For obtaining a 
most nearly geometric sequence of approximations, assume that 

U(x.e) ~ U|"(x) + cLIt (x) + e^U 3 (x) + . . . . (2.6) 

Successive approximations to U(x) are then defined by n-term truncations of the series (2.6); 

n 

Un = ^ c'-' Ui'(x) (2.7) 

i=I 

Although (2.6) is equivalent to a Taylor series or asymptotic series expansion about e = 0, its conver- 
gence or lack of convergence at e = 1 is not of particular significance for applicability of eq. (2.3) (see 
[3] ). If the series (2.6) is substituted into the problem of eq. (2.4) and condition (2.2) and coefficients 

t 

of powers of e are collected, one obtains equations to solve for the Un : 


LUj'-F - NUo in R; BU]' = G(x) on B ; (2.8a) 

LU2” = NUi'-NUo in R; BLI2' - 0 on B ; (2.8b) 

LU3' = N2' lUn'X'i'J iiiR; BU3' - 0 on B; (2.8c) 

i 

in which N2 is defined by the perturbation expansion 

NU = NUi' + eN2' ]U2 ’Ui'f +0(e-). (2.9) 

With the definitions (2.7) and 

N2 ]U2>Ui} = NU]' + eN2' ]U2',Ui'[ (2.10) 

one can also solve the following equations for the successive approximations, 1)^; 

LUi -F = NUo R: BUi = G(x) on B (2.11a) 

LU2-F - NUi in R; BU2 - G(x) on B (2.11b) 

LU3-F = N2 jUo.Uif in R; BU3 = Ci(x) on B (2.11c) 


in which it has been assumed that e- 1. Note that if the right side of eq. (2.1) is linear in U(x), then 
the problems for the successive in eqs. (2.11) are the same as (2.5) for Method 1 . 

We denote as Method 2(a) the solution of eqs. (2.1 1) for three successive iterates and substitution 
of the results for one component of each 1)^ into (2.3) to obtain an improved approximation. (If NU 
is linear in U. this is the same as Method 1(b)). Note that when the solution is near to convergence at 
any x, significant errors will be introduced by the loss of significant figures in applying eq. (2.3). ' 

An alternative procedure (denoted as Method 2(b)) that eliminates tiie difficulty near convergence 
is to replace eq. (2.3) by the equivalent expression (at e = 1 ): 

, (U2 )“ 

u*(x) = U] - — r , (2.12) 

113 - U2 

where each u^ (x) is a single component of the vector Up (x). That is, eqs. (2.8) are solved for U^(x), 
and (2.12) is used for extrapolation. 

In a numerical solution, u*(x) can be used as the next Uq(x) in a repetition of the sequence. 

3. Example Problems and Comparison of Metliods 

This section gives simple analytical one-dimensional examples for illustration and comparison of 
the methods. 

3.1 Example 1. Consider the nonlinear problem 

(d/dx + 1)u = (l/2)u“ inO<x<<», (3.1a) 

u(0 ) = 1 . (3.1b) 
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The iterative solution by Method 1 is found from 


(d/dx+ i)Un = > MO) = I • (3.2) 

The analytical solutions for n= 1.2,3 (assuming Uq = 0) are; 



U|(x) = e'-’^ , 

(3.3a) 


UTf X) = e'^ [ 1 + p( x)J . 

(3.3b) 


U 3 (x) - e'-’^ ( 1 + p(x) + p^(.\) + j P'^(x)| . 

(3.3c) 

where 


p(x) = ( 1 / 2 ) ( 1 -e*^). 

(3.4) 

For Method 2, 

the artificially extended equation is; 



(d/dx+l)u = ( 1 -e) (l/ 2 )Uo^ + e(l/ 2 )u“ , 

(3.5 a) 


u( 0 ) - 1 . 

(3.5bX 

Substitution of 


$ ! 'J * 

u = u j (x) + 6 U 2 (x) + e^U 3 (x) + . . . 

(3.6) 

into (3.5) leads to 


(d/dx+ l)ui' = (l/2)uo^. uj'(O) = 1 , 

(3.7a) 


(d/dx + Din = ( 1 / 2 ) |(u] )“-Uo^l , in ( 0 ) = 0 , 

(3.7b) 


(d/dx + 1 )U 3 = U| UT , 1 ^ 3 ( 0 ) ~ 0 > 

(3.7c) 

or equivalently, with e = 1 and eq. (2.7), 



(d/dx + Dll] = (l/ 2 )Uo-. ui( 0 ) = 1 , 

(3.8a) 


(d/dx+l)u 2 == (l/ 2 )uj“, infO) = 1 , 

(3.8b) 


(d/dx + l)u 3 = ( 1 / 2 ) 02 " -(l/ 2 )(ui-U 2 )-, 113 ( 0 ) “ * • 

(3.8c) 

The analytical solutions to (3.7) with Uq = 0 are: 



Un’(x) = e'’‘[p(x)]"'^ , 

(3.9) 


where p(x) is given by (3.4) and where the solutions Un to (3.8) are given by (2.7). Evaluations of 
these solutions at x = 1 and applications of the appropriate forms of the Aitken/Shanks formula are 
given in Table 1. The results for the extrapolated solution u* may be compared with the exact solution 



Table 1. Results of Example 1 at x = 1 (Uq = 0) 


METHOD; 

Kb) 

2(a) 

2(b) 

EQUATIONS: 

(3.2) & (2.3) 

(3.8) & (2.3) 

(3.7) & (2.12) 

n 

Un(l) 

Un(l) 


1 

0.3678794412 

0.3678794412 

0.3678794412 

■y 

.4841515202 

.4841515202 

.2325441579 

3 

.5247721376 

.5209005060 

.1469959430 

u*(l) = 

.546583145 

.537882842 

.5378828426 

Exact u(l) = 

.5378828428 

.5378828428 

.5378828428 


evaluated at x=l: u(l) = 0.5378828428 to ten significant figures. Wc note first that the extrapolated 
solution u* by Method 1(b) is somewhat closer to the exact value than 113. but not significantly closer. 
We note further that the third approximation, U3, by Method 2(a) is not as good an approximation as 
U3 in Method 1, but that the extrapolated solutions by Methods 2(a) and 2(b) are exact except for loss 
of 1 or 2 significant figures. (Method 2(a) is less exact because of loss of significant figures in (2.3).) 

The striking accuracy of Method 2 in this example occurs because the sequence of solutions produced by 
Method 2 is precisely geometric, i.e. Un+i/u^ = constant for all n at a given x. The difference from 
Method 1 is seen by comparing cqs. (3.2) with (3.8), in which (3.8c) has an additional term that produces 
the geometric sequence. 

3.2 Example 2. Consider next an example which is linear (so that .Method 2(a) would give the same 
results as Method 1(b)), but for which the iterative sequence is “nearly geometric.” Let us use Method 
2(b) for this example (eqs. (2.8) with (2.6), (2.7), and (2.12)). 

The problem is 


-^-2=-x^ -2u in 0<x<<» . u(0) = 0 , (3.11) 

dx dx 

which is written in this way in analogy to more complex problems in which one may put a very simple 
operator on the left and the rest of tlie terms on the right for iteration. (One can also shift the term 2u 
to the left side, with very similar results.) Tiie artificially extended equation is 

du fill 

^ - 2 = (5-e)(-x-^ -2uo) + e(-x^ - 2u) in 0<x<«> 

u(0) = 0 


(3.12) 
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Substitution of (3.6) leads to (with Uq - 0): 


du]'/dx-2 = 0. ui'(0) = 0. (3.13a) 

din’/dx - -xdui /dx- 2 uj , ut (0) = 0, (3.13b) 

duj ,'dx = -X duj) /dx - 2 u 2 , uj (0) = 0 . (3.13c) 

The analytical solutions are 

Un’(x) = (-l)”"^' (n+l)x” “ (3.14) 


and the successive approximations are given by (2.7). The sequence (3.14) is not geometric, but since 

/ t 

2imn-K» lun+](x)/un (x)l exists at given x. the sequence is “nearly geometric” (3] . Evaluation of the 

It* 

solutions (3.14) at x = 0.5 gives (uj , ut , U 3 ) = (1 .00. -.75, .50) so that the successive approximations 

f 

are (uj, 03 , 03 ) = (1.00. .25, .75). Substitution of the Up into (2.12) gives u*(.5) = 0.55, which com- 
pares well with the exact solution to (3.1 1), 


u(x) = 


(2x + X-) ( 1 + x)‘^ , 


(3.15) 


frofti which u(0.5) = 5/9 0.555555 .... 

4. Transonic Flow Over an Airfoil 

For application of the methods described above, consider two-dimensional, steady, inviscid flow 
over a thin symmetrical parabolic-arc airfoil in a subsonic free stream. At high subsonic Mach numbers, 
part of the flow can be supersonic, so we consider the transonic small-disturbance equations, which are 
nonlinear elliptic partial differential equations in subsonic regions and hyperbolic equations in super- 
sonic regions. Transition of the velocity field from a subsonic region to the embedded supersonic zone 
is smooth, but transition from the supersonic to subsonic region is usually discontinuous, through a 
shock wave. The improved finite-difference method of Murman and Cole [ 19-22) captures the shock 
waves (in a fully conservative way) but spreads the rapid transition over several mesh points. 

In [ 8 ] a semidirect finite-difference method, based on the use of a fast direct Cauchy-Riemann 
solver [15], was applied to solving the equivalent of Murman ’s transonic finite-difference equations [21 ] 
iteratively for the perturbation velocities, u and v. (The iteration procedure has been formulated in 
such a way that at nonelliptic points terms on the right side of the difference equations cancel out the 
elliptic character of the left side when the iterated solution converges.) Both Methods 1(a) and 2(a) 
described above worked well lor subcritical and for slightly supercritical (local Mach number > 1 ) flows, 
except that Method 2(a) could be used only before any part of the solution was nearly converged. In 
[9] the method was extended to strongly supercritical flow by the addition of stabilizing terms to the 
difference equations and lo the (’auchy-Rieinann solver | 18|. Also iniroiluced in [9| was (he method 
version denoted here as Method 2(h), which can he used when the solution is nearly converged. In 
smooth subsonic flows the acceleration technique is effectively used repeatedly. However, in transonic 
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flows with strong shock waves, the acceleration technique is not helpful at the beginning of the iteration 
when the shock wave and its location are not well defined. Therefore in (9) it was considered desirable 
to use the straightforward iteration Method 1(a) until the maximum residual is reasonably small, so that 
the supersonic region is nearly defined, and then use Method 2(b) to extrapolate three iterates to a final 
solution. A fully conservative second-order-accurate formulation has been introduced in [ 10] , and so a 
fomulation that includes either Murman’s fully-conservative first-order-accurate formulation or the 
second-order formulation will be used here. 

4. 1 Governing equations and boundary conditions. Let the dimensionless X and Y axes be 
respectively along and normal to the airfoil chord, the free-stream Mach number be M^q < 1 , and the 
dimensionless velocity components in the X and Y directions be U,V. One may then define perturba- 
tion velocity components u.v through a Prandtl-Glauert transformation with (3 = (1 - 

U = 1 +(t//3)u, V = rv, Y = y//3 , X = x, (4.1) 

which amounts to shifting and scaling of certain terms (cf 1 17, 8-10] , so that the transonic small 
disturbance equations take the form 

&y ~ 0 • i^y ~ 'x “ 0 (4.2a,b) 

where 

f =■ flu) = u - au^ . g = g(v) = V, (4.3a,b) 

a = r(7+l)M^/2()3, (4.4) 

in which a is a transonic similarity parameter and r is an airfoil thickness ratio. Eqs. (4.2) are often 
written in terms of a perturbation velocity potential <j> defined by u = ~ 0y> and all the develop- 

ments to be described apply as well to that potential equation. 

The equation system (4.2) is elliptic, parabolic, or hyperbolic depending on whether u - is 
negative, zero, or positive, where the transformed critical velocity is u^-j^ = l/2a. The corresponding 
pressure coefficient is Cp = -2 (t/(3)u. 

The linearized surface boundary condition for the symmetrical parabolic-arc airfoil, whose upper 
surface is given by Y 5 (x) = rF(x) = r(0.5 - 2x") in -.5 < x < .5 (with F(x) = 0 in |x| > .5 ), and the 
conditions at infinity are 

v(x,0+) = f'(x) , (4.5a) 

u,v-*-0 as x^ + y (4.5b) 

Eq. (4.2a) is written in a “conservation-law” (or divergence) form, in terms of flux components f 
and g. Therefore discretized forms of (4.2a), for numerical solution, can represent in a fully conserva- 
tive way either that differential equation or the corresponding integral form. These discretized forms 
can thus be formulated correctly to represent transitions between elliptic and hyperbolic regions [21,10] . 
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Since only the term in the system (4.2) determines the type of point (depending on the local 
value of u), one can write the genera! type-dependent difference equations in the form: 

(fx)T + (gy)c = 0. (i>y)C ” ~ 0- (4.6a,b) 

where subscript C indicates a central-differenced representation of a derivative and subscript T, which 
indicates type-dependent differencing, may be replaced by E, H. P, or S at points defined respectively 
as elliptic, hyperbolic, parabolic, or shock points [21,10] . At all points where the difference equations 
are clearly elliptic or hyperbolic, subscripts E or H are used. Transition points from elliptic to hyper- 
bolic (progressing downstream from left to right) are P points, and transitions from hyperbolic to 
elliptic are S points. 

For defining the finite-difference operators. Fig. 1 shows a staggered u,v mesh, with the shaded 
area indicating a mesh cell for eq. (4.2a). The center of a mesh cell is the point at which <j> would be 
defined on a conventional mesh and is the point that is designated E, FI, P, or S. The indices j and k 
indicate respectively the x and y directions. Second-order-accurate centra! differences are 

(“x)c = - ''j-1 • (''y)c = (''j.k " Vj.k-1 )/^y - ) 

(Uy)(; = (iij k-H - Uj,k)/Ay , = (vj+i k-vj k)/Ax. ) 

In general, (f^)j is represented by 

Ax(fx)T = Alj,k = "(fcVl k (4.8a) 

where 

and where uq is either a “hyperbolic form” uj^ or an “elliptic form” U£. With (i) the definition (4.8) 
for the difference operator (fx) 7 - G‘) a condition to determine whether each (u(^)j is represented by 
U£ or ,ujj, and (iii) specifications of uj? and upj to obtain the finite differences (4.8a) to the order of 



I AX 


Fig. 1 — Differencing mesh and mesh cell. 
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accuracy desired, all four type-dependent operators are obtained. As derived in 1 10) , (i) the second- 
order-accuratc elliptic operator, (ii) either the tlrst-order or seeond-order-aceurate hyperbolic operator, 
and (iii) the correspondinu parabolic-point and shock-point operators are all produced in (4.6) with (4.7) 
and (4.8) by the following relatiotiship; 

'I'G^.i.k ^ <'-''j.k»»j.k-^®j.k(^“j-l.k + <'-^)“j-2.kl ’ 

where 

Oj k = 0 (and UG = U£) if U|^ < Uj;.^ , 

= 1 (and uq = Uj,{) if uj > u^j^ , 

"j.k ^ ^“j,k+^"j-l,k>/<'+^>- 

and where X=1 for the tlrst-order-accurate hyperbolic operator, X = 2 for the second-order-accurate 
hyperbolic operator, and 5 is a parameter that may be varied from 0 to o° but is derived as unity for 
Murman’s first-order-accurate operators (21 ). As an example to illustrate, suppose X = 1, 5 = 1, and 
Uj k < 2nd Uj.| k > Then for the shaded mesh cell in Fig. 1 , eqs. (4.8) - - (4. 10) give 

^^j,k ~ “j,k ~ *^j-2,k “ ~ *^f-2.k^ 

which is equivalent to Murman's (2 1 1 first-order shock-point operator. In a similar way the Krupp- 
Murman first-order parabolic operator )20) is also obtained. Both the first- and second-order-accurate 
hyperbolic operators given by (4.8) - (4. 10) with X = 1 and X = 2 are equivalent to upwind difference 
operators originally propo.sod by Murman and Cole ) 19] ; the fully conservative second-order P and S 

operators were introduced in ( 10] . Analysis of all these E, H, P, and S operators [ 1 0] has verified 

their consistency, accuracy, and stability in the examples computed. 

Because of the slow iterative convergence of the second-order-accurate iterative method to be 
described, two methods of adding artificial viscosity have been proposed and used [ 10] . Both leave 
the scheme fully conservative and formally second-order-accurate. 

The boundary conditions for the finite-difference equations (4.6) are the same as (4.5) but with 
(4.5b) replaced by a far-field condition on an outer rectangular boundaiy B; 

Uj k = ug(.\.y) or vj k = vg(x,y) on B (4.12) 

where, for example, ug and vg are given by a Prandtl-Glauert solution (see [ 15,8,9] ). 

For solution of eqs. (4.6) with (4.7) through (4.1 1) and with conditions (4.12) by the semidirect 
methods, one must rearrange the equations so that the left side is an appropriate elliptic operator and 
provides a stable iteration scheme. One first adds (Ux^C~^^x^T to both sides of (4.6a) to obtain • 

^“xk: ^''y^C “ ^^xV’ “ ^^x^T ’ (4.13a) 

(Uy)c-(vx)c = 0. (4.13b) 


(4.10) 


1 n 
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This set contains a central-dilTerenced elliptic operator on the left side regardless of the local type of the 
equations. The nonlinear type-dependent term has been shifted to the right side where, in an iterative 
procedure, it can be computed from a previous iteration. Although the iteration of these equations (8) 
converged well for subsonic and slightly supercritical How, it was found (9,10) that terms with param- 
eters multiplying Uj and Uj_j j. needed to be added to both sides of (4.13a) to produce iterative 
convergence at higher Mach numbers. more specific form of the difference equations, in which the 
second-order-accurate relations ( 4 . 7 ) have been substituted, is 


Dj k(ihv) = (4.14a,b) 

in which 

Dj kfii-v) = (1 -a|)iij,k “(1 +M*'(vj,k - Vj k-l) . (4.15a) 

Ej,k(ihv) H (Uj k+i-uj k)-Mvj+],k-''j,k>> <415b) 

Rj.k(u) = f‘-“l>“j.k-('+“2)“j-l,k-^fj.k • (4.15c) 

and where Atj k is defined by ecjs. (4.8) (4. 1 1 ) and p = Ay/Ax. The formal order of accuracy of 

eqs. (4.14) depends on the value of .\ used in (4.9). 

4.2 Equations for Method 1(a). As described in section 2.2 above, the straiglitforward iteration 
Method 1(a) for eqs. (4. 14) is simply 

Dj.k(^n-'^n) = •^j,k<»n-l >’ kfu^.v^) = 0 . (4.16a,b) 

For determining each oj k in (4.10), eq. (4.1 1) uses Uj^. j. The presence of ajUj k and a-jUj.j k on 
both sides of eq. (4.16a) allows the interpretation and treatment of these terms as an off-centered time 
derivative, 3u/9t, multiplied by a constant. When the solution converges, these terms cancel out. The 
semidirect Method 1(a) proceeds by solving the left side of (4.16) in terms of the known right side by an 
“extended Cauchy-Riemann” solver [ 181 for Uj^ and Vj^ at all points simultaneously. The iteration with 
aj or U 2^0 needs a reasonable (but very roughly approximate) initial approximation (Uq), such as a 
Prandtl-Glauert solution. Ref. [lO] gives variable specifications of a-> for best convergence. 


The boundary conditions on (4.1 6) arc 

vn(x,0+) = F'(x) . (4.17a) 

“n ~ *-*B ''n ~ ''B on B . (4.17b) 

4.3 Equations for Method 2(b). The artificially extended form. (2.4), of eqs. (4.14) is 

Dj k(u,v) = (1 -e)Rj k(U(,) + eRj k(u) , (4.18a) 

Ejk(u.v) = 0. (4.18b) 
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For Method 2(b) assume that 


u(x,y,e) = u/(x,y) + eu7Vx,y) + e^UoXx.y) + ... , (4.19a) 

v(x.y.e) =vi'(x,y) + ev','(x.y) + e"V3(x,y) + ... . (4.19b) 

The successive approximations are then (lor n = 1,2,3 . . .) 


U|^ -- ^ Uj\'x,y) . Vp - ^ ^ e'''vj'(x,y) . (4.20) 

i=l i-l 

Substitution of (4.19) into (4.18) leads to 

= *^n-l • Ejk(u„',v„) = 0, (4.21) 

where: 

Ro = Rj,k(u„) (4.22a) 

Rj = Rj|^(u| )~ Rj k^*^o^ (4.22b) 

(with Uq being used in (4.1 1) in determining oj for use in Rj i^luj )) and 

R 2 = (1 -ai)(u 2 ')j,k +“2>^^'2’^j,k '(^‘2)j,k ’ (4.22c) 

Af2 = (f2)j.k-(‘2>M,k' 

(f2)j^k = - <^j,k) (^»:'>j,k - -=»<«r“2'>j,k J 

+ Oj i^ I [Mii2')j-i,k 

- 2a(\(Uj')j.j |. + (l-X)(U]')j_ 2 Vl .k ^ * ~^^*^'2'V2,k^f • (4.23b) 

The boundar>' conditions are: 

viVx.O"^) = FXx) ; Vj.^'(x,0"'’) = 0 (n = 2,3) : (4.24a) 

Uj' = Ur or ' ) = vr on B; (4.24b) 

Up =0 or Vp' = 0 on B (n = 2,3) . (4.24c) 


With some reasonable approximation for (Uq)j j.. such as a nearly converged solution by Method 1 (a), 
eqs. (4.21), with n= 1.2,3, give three successive approximations U| , 02 , 113 ’ at each j,k to use in 
( 2.1 2 ) to obtain an extrapolated solution. 

4.4 Results and discussion. A research computer program written to solve the transonic small 
disturbance equations by the methods described above for a biconvex airfoil at zero incidence, includes 
the option of switching after some iterations by Method 1(a) to the extrapolation technique. Method 
2(a). A conversational version of the program, for interacting with the program, was run on an IBM 
360/67 computer, and computing times were measured on a Control Data 7600 computer. 
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Pressure distributions have been computed for a range of subsonic and transonic Mach numbers 
from both first- and second-order-accuratc formulations. Examples by Method 1(a) are shown on Fig. 2 
for a thickness ratio of 1 0 percent and M^^ = 0.825. For this calculation the boundaries were at one-half 
chord upstream and downstream of the airfoil edges and at 3.5 chords above the airfoil. The results com- 
puted on a 39X32 uniform mesh compare well with a line-relaxation program (22) , which uses a variable 
and finer mesh. On a very' coarse ( 19X32) mesh, with only 10 mesh intervals on the airfoil chord, the 
first-order-accurate results, of course, are not good. The shock is badly smeared, and an anomalous jump 
behind the sonic point that is characteristic of the first-order P operator is exaggerated on the coarse 
mesh. However, the second-order-accurate results arc very smooth through the sonic point and are 
surprisingly accurate. 




Moo = 0.825. r=O.IO 

o 1st ORDER. 
19x32 MESH 

A 2nd ORDER. 
19x32 MESH 


MURMAN et Ql, 

VARIABLE MESH. 
REF. 22 


Fig. 2 — Pressure on a thin biconvex airfoil. 
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Figure 3 shows an interesting el feet oT switching to Method 2 bet'ore the iteration has converged 
enough, when the types ol all points are not yet t|uite the same as the final types. Method 1(a) was used 
for nine iterations; then Metlu)d 2(b) was used to obtain the three successive terms at each point and the 
extrapolated solution shown in Fig. 3. A property of the Aitken /Shanks extrapolation as used in 
Method 2 is that all the significant figures of tlie three successive approximations at any point contain 
information about the exact solution, even though those successive approximations themselves are not 
very close to the exact solution (see example problems above, in section 3). It thus appears possible in 
Fig. 3 that this procedure may be picking up the fact that the exact solution to the equations (or the 
solution on a very fine mesh) has the well-known logarithmic singularity just behind the shock, even 
though the. converged solution on the coarse mesh smears over this singularity. Even the finer mesh used 
by tlie program in (221 was not fine enough to pick up the singularity, partly because that point appar- 
ently occurs between the mesh points for this case. This phenomenon illustrated in Fig. 3 is not an iso- 
lated case but is a typical occurrence in Method 2. It may be that the numerical solution in Fig. 3 is as 
good as representation of the exact solution to the equations as is the fully converged solution in Fig. 2(a) 
(circles). 

The most significant property of the semidirect method is the relatively short computing time 
required. On the 39X32 mesh, the time per iteration was measured as 40 milliseconds in a very ineffici- 
ently coded program, but for various reasons discussed in [10] it is expected to be reduced to 20 ms. 

(Tlie direct solver requires only 14 ms) The subcritical cases were sufficiently converged in 3 iterations 
or less, and a slightly supercritical case (first-.order-accurate, using Method 2) required 6 iterations. The 



X 


Fig. 3 — Pressure distribution resulting from Aitken/Shanks extrapolation (Method 2(b)) 
before iterative convergence. 


s,._ 



first-order-accurate case shown in Fig. 2(a) required 20 iterations by Method 1 and, as described above, 
tlie results of Fig. 3 required only 9 iterations by Method 1(a) followed by 3 more by Method 2(b). 

At this writing, tlie program has not yet been written for the above formulation that includes the 
second-order-accurate fonnulation in Method 2. It is e.xpected that when this is done, the program can 
be run rapidly with the first-order (X = 1) operators on the very coarse mesh using Method 1, then 
switched to second-order (X = 2) and Method 2 for final e.xtrapolation. 


5. Concluding Remarks 

It has been shown that a special procedure (Method 2) is effective for obtaining most appropriate 
successive appro.ximations for use in the Aitken extrapolation formula for accelerating the iterative con- 
vergence of numerical solutions to nonlinear partial differential equations. The procedure is based on the 
combined use of artificial perturbation-series expansions and an artificially extended equation. It was 
shown in a previous paper [8] that one version of the technique was very effective for accelerating itera- 
tive convergence when the solutions are smooth. The method, in a modified version, has now been applied 

with some success to a strongly supercritical transonic flow problem, in which the flow equations are of 

♦ 

mixed type and whose solutions have shock-wave discontinuities. The method is expected to be extended 
to more general flows, including lifting airfoils and three-dimensional flows. 
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